Aminet 1 (Walnut Creek)

home *** CD-ROM | disk | FTP | other *** search

/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / usenet / sources / volume2 / aplictns / matlab / src.6 < prev next >

Wrap

Internet Message Format | 1988-11-02 | 54KB

Path: xanth!nic.MR.NET!hal!cwjcc!mailrus!bbn!ulowell!page From: page@swan.ulowell.edu (Bob Page) Newsgroups: comp.sources.amiga Subject: v02i046: matlab - matrix laboratory, Part06/11 Message-ID: <10021@swan.ulowell.edu> Date: 2 Nov 88 21:45:11 GMT Organization: University of Lowell, Computer Science Dept. Lines: 1220 Approved: page@swan.ulowell.edu Submitted-by: strovink%galaxy-43@afit-ab.arpa (Mark A. Strovink) Posting-number: Volume 2, Issue 46 Archive-name: applications/matlab/src.6 # This is a shell archive. # Remove everything above and including the cut line. # Then run the rest of the file through sh. #----cut here-----cut here-----cut here-----cut here----# #!/bin/sh # shar: Shell Archiver # Run the following text with /bin/sh to create: # src-6 # This archive created: Wed Nov 2 16:20:42 1988 cat << \SHAR_EOF > src-6 SMALL = 1.D0/2.D0**48 C C DETERMINE WHAT IS TO BE COMPUTED. C WANTU = .FALSE. WANTV = .FALSE. JOBU = MOD(JOB,100)/10 NCU = N IF (JOBU .GT. 1) NCU = MIN0(N,P) IF (JOBU .NE. 0) WANTU = .TRUE. IF (MOD(JOB,10) .NE. 0) WANTV = .TRUE. C C REDUCE X TO BIDIAGONAL FORM, STORING THE DIAGONAL ELEMENTS C IN S AND THE SUPER-DIAGONAL ELEMENTS IN E. C INFO = 0 NCT = MIN0(N-1,P) NRT = MAX0(0,MIN0(P-2,N)) LU = MAX0(NCT,NRT) IF (LU .LT. 1) GO TO 190 DO 180 L = 1, LU LP1 = L + 1 IF (L .GT. NCT) GO TO 30 C C COMPUTE THE TRANSFORMATION FOR THE L-TH COLUMN AND C PLACE THE L-TH DIAGONAL IN S(L). C SR(L) = WNRM2(N-L+1,XR(L,L),XI(L,L),1) SI(L) = 0.0D0 IF (CABS1(SR(L),SI(L)) .EQ. 0.0D0) GO TO 20 IF (CABS1(XR(L,L),XI(L,L)) .EQ. 0.0D0) GO TO 10 CALL WSIGN(SR(L),SI(L),XR(L,L),XI(L,L),SR(L),SI(L)) 10 CONTINUE CALL WDIV(1.0D0,0.0D0,SR(L),SI(L),TR,TI) CALL WSCAL(N-L+1,TR,TI,XR(L,L),XI(L,L),1) XR(L,L) = FLOP(1.0D0 + XR(L,L)) 20 CONTINUE SR(L) = -SR(L) SI(L) = -SI(L) 30 CONTINUE IF (P .LT. LP1) GO TO 60 DO 50 J = LP1, P IF (L .GT. NCT) GO TO 40 IF (CABS1(SR(L),SI(L)) .EQ. 0.0D0) GO TO 40 C C APPLY THE TRANSFORMATION. C TR = -WDOTCR(N-L+1,XR(L,L),XI(L,L),1,XR(L,J),XI(L,J),1) TI = -WDOTCI(N-L+1,XR(L,L),XI(L,L),1,XR(L,J),XI(L,J),1) CALL WDIV(TR,TI,XR(L,L),XI(L,L),TR,TI) CALL WAXPY(N-L+1,TR,TI,XR(L,L),XI(L,L),1,XR(L,J), * XI(L,J),1) 40 CONTINUE C C PLACE THE L-TH ROW OF X INTO E FOR THE C SUBSEQUENT CALCULATION OF THE ROW TRANSFORMATION. C ER(J) = XR(L,J) EI(J) = -XI(L,J) 50 CONTINUE 60 CONTINUE IF (.NOT.WANTU .OR. L .GT. NCT) GO TO 80 C C PLACE THE TRANSFORMATION IN U FOR SUBSEQUENT BACK C MULTIPLICATION. C DO 70 I = L, N UR(I,L) = XR(I,L) UI(I,L) = XI(I,L) 70 CONTINUE 80 CONTINUE IF (L .GT. NRT) GO TO 170 C C COMPUTE THE L-TH ROW TRANSFORMATION AND PLACE THE C L-TH SUPER-DIAGONAL IN E(L). C ER(L) = WNRM2(P-L,ER(LP1),EI(LP1),1) EI(L) = 0.0D0 IF (CABS1(ER(L),EI(L)) .EQ. 0.0D0) GO TO 100 IF (CABS1(ER(LP1),EI(LP1)) .EQ. 0.0D0) GO TO 90 CALL WSIGN(ER(L),EI(L),ER(LP1),EI(LP1),ER(L),EI(L)) 90 CONTINUE CALL WDIV(1.0D0,0.0D0,ER(L),EI(L),TR,TI) CALL WSCAL(P-L,TR,TI,ER(LP1),EI(LP1),1) ER(LP1) = FLOP(1.0D0 + ER(LP1)) 100 CONTINUE ER(L) = -ER(L) EI(L) = +EI(L) IF (LP1 .GT. N .OR. CABS1(ER(L),EI(L)) .EQ. 0.0D0) * GO TO 140 C C APPLY THE TRANSFORMATION. C DO 110 I = LP1, N WORKR(I) = 0.0D0 WORKI(I) = 0.0D0 110 CONTINUE DO 120 J = LP1, P CALL WAXPY(N-L,ER(J),EI(J),XR(LP1,J),XI(LP1,J),1, * WORKR(LP1),WORKI(LP1),1) 120 CONTINUE DO 130 J = LP1, P CALL WDIV(-ER(J),-EI(J),ER(LP1),EI(LP1),TR,TI) CALL WAXPY(N-L,TR,-TI,WORKR(LP1),WORKI(LP1),1, * XR(LP1,J),XI(LP1,J),1) 130 CONTINUE 140 CONTINUE IF (.NOT.WANTV) GO TO 160 C C PLACE THE TRANSFORMATION IN V FOR SUBSEQUENT C BACK MULTIPLICATION. C DO 150 I = LP1, P VR(I,L) = ER(I) VI(I,L) = EI(I) 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE 190 CONTINUE C C SET UP THE FINAL BIDIAGONAL MATRIX OR ORDER M. C M = MIN0(P,N+1) NCTP1 = NCT + 1 NRTP1 = NRT + 1 IF (NCT .GE. P) GO TO 200 SR(NCTP1) = XR(NCTP1,NCTP1) SI(NCTP1) = XI(NCTP1,NCTP1) 200 CONTINUE IF (N .GE. M) GO TO 210 SR(M) = 0.0D0 SI(M) = 0.0D0 210 CONTINUE IF (NRTP1 .GE. M) GO TO 220 ER(NRTP1) = XR(NRTP1,M) EI(NRTP1) = XI(NRTP1,M) 220 CONTINUE ER(M) = 0.0D0 EI(M) = 0.0D0 C C IF REQUIRED, GENERATE U. C IF (.NOT.WANTU) GO TO 350 IF (NCU .LT. NCTP1) GO TO 250 DO 240 J = NCTP1, NCU DO 230 I = 1, N UR(I,J) = 0.0D0 UI(I,J) = 0.0D0 230 CONTINUE UR(J,J) = 1.0D0 UI(J,J) = 0.0D0 240 CONTINUE 250 CONTINUE IF (NCT .LT. 1) GO TO 340 DO 330 LL = 1, NCT L = NCT - LL + 1 IF (CABS1(SR(L),SI(L)) .EQ. 0.0D0) GO TO 300 LP1 = L + 1 IF (NCU .LT. LP1) GO TO 270 DO 260 J = LP1, NCU TR = -WDOTCR(N-L+1,UR(L,L),UI(L,L),1,UR(L,J), * UI(L,J),1) TI = -WDOTCI(N-L+1,UR(L,L),UI(L,L),1,UR(L,J), * UI(L,J),1) CALL WDIV(TR,TI,UR(L,L),UI(L,L),TR,TI) CALL WAXPY(N-L+1,TR,TI,UR(L,L),UI(L,L),1,UR(L,J), * UI(L,J),1) 260 CONTINUE 270 CONTINUE CALL WRSCAL(N-L+1,-1.0D0,UR(L,L),UI(L,L),1) UR(L,L) = FLOP(1.0D0 + UR(L,L)) LM1 = L - 1 IF (LM1 .LT. 1) GO TO 290 DO 280 I = 1, LM1 UR(I,L) = 0.0D0 UI(I,L) = 0.0D0 280 CONTINUE 290 CONTINUE GO TO 320 300 CONTINUE DO 310 I = 1, N UR(I,L) = 0.0D0 UI(I,L) = 0.0D0 310 CONTINUE UR(L,L) = 1.0D0 UI(L,L) = 0.0D0 320 CONTINUE 330 CONTINUE 340 CONTINUE 350 CONTINUE C C IF IT IS REQUIRED, GENERATE V. C IF (.NOT.WANTV) GO TO 400 DO 390 LL = 1, P L = P - LL + 1 LP1 = L + 1 IF (L .GT. NRT) GO TO 370 IF (CABS1(ER(L),EI(L)) .EQ. 0.0D0) GO TO 370 DO 360 J = LP1, P TR = -WDOTCR(P-L,VR(LP1,L),VI(LP1,L),1,VR(LP1,J), * VI(LP1,J),1) TI = -WDOTCI(P-L,VR(LP1,L),VI(LP1,L),1,VR(LP1,J), * VI(LP1,J),1) CALL WDIV(TR,TI,VR(LP1,L),VI(LP1,L),TR,TI) CALL WAXPY(P-L,TR,TI,VR(LP1,L),VI(LP1,L),1,VR(LP1,J), * VI(LP1,J),1) 360 CONTINUE 370 CONTINUE DO 380 I = 1, P VR(I,L) = 0.0D0 VI(I,L) = 0.0D0 380 CONTINUE VR(L,L) = 1.0D0 VI(L,L) = 0.0D0 390 CONTINUE 400 CONTINUE C C TRANSFORM S AND E SO THAT THEY ARE REAL. C DO 420 I = 1, M TR = PYTHAG(SR(I),SI(I)) IF (TR .EQ. 0.0D0) GO TO 405 RR = SR(I)/TR RI = SI(I)/TR SR(I) = TR SI(I) = 0.0D0 IF (I .LT. M) CALL WDIV(ER(I),EI(I),RR,RI,ER(I),EI(I)) IF (WANTU) CALL WSCAL(N,RR,RI,UR(1,I),UI(1,I),1) 405 CONTINUE C ...EXIT IF (I .EQ. M) GO TO 430 TR = PYTHAG(ER(I),EI(I)) IF (TR .EQ. 0.0D0) GO TO 410 CALL WDIV(TR,0.0D0,ER(I),EI(I),RR,RI) ER(I) = TR EI(I) = 0.0D0 CALL WMUL(SR(I+1),SI(I+1),RR,RI,SR(I+1),SI(I+1)) IF (WANTV) CALL WSCAL(P,RR,RI,VR(1,I+1),VI(1,I+1),1) 410 CONTINUE 420 CONTINUE 430 CONTINUE C C MAIN ITERATION LOOP FOR THE SINGULAR VALUES. C MM = M ITER = 0 440 CONTINUE C C QUIT IF ALL THE SINGULAR VALUES HAVE BEEN FOUND. C C ...EXIT IF (M .EQ. 0) GO TO 700 C C IF TOO MANY ITERATIONS HAVE BEEN PERFORMED, SET C FLAG AND RETURN. C IF (ITER .LT. MAXIT) GO TO 450 INFO = M C ......EXIT GO TO 700 450 CONTINUE C C THIS SECTION OF THE PROGRAM INSPECTS FOR C NEGLIGIBLE ELEMENTS IN THE S AND E ARRAYS. ON C COMPLETION THE VARIABLE KASE IS SET AS FOLLOWS. C C KASE = 1 IF SR(M) AND ER(L-1) ARE NEGLIGIBLE AND L.LT.M C KASE = 2 IF SR(L) IS NEGLIGIBLE AND L.LT.M C KASE = 3 IF ER(L-1) IS NEGLIGIBLE, L.LT.M, AND C SR(L), ..., SR(M) ARE NOT NEGLIGIBLE (QR STEP). C KASE = 4 IF ER(M-1) IS NEGLIGIBLE (CONVERGENCE). C DO 470 LL = 1, M L = M - LL C ...EXIT IF (L .EQ. 0) GO TO 480 TEST = FLOP(DABS(SR(L)) + DABS(SR(L+1))) ZTEST = FLOP(TEST + DABS(ER(L))/2.0D0) IF (SMALL*ZTEST .NE. SMALL*TEST) GO TO 460 ER(L) = 0.0D0 C ......EXIT GO TO 480 460 CONTINUE 470 CONTINUE 480 CONTINUE IF (L .NE. M - 1) GO TO 490 KASE = 4 GO TO 560 490 CONTINUE LP1 = L + 1 MP1 = M + 1 DO 510 LLS = LP1, MP1 LS = M - LLS + LP1 C ...EXIT IF (LS .EQ. L) GO TO 520 TEST = 0.0D0 IF (LS .NE. M) TEST = FLOP(TEST + DABS(ER(LS))) IF (LS .NE. L + 1) TEST = FLOP(TEST + DABS(ER(LS-1))) ZTEST = FLOP(TEST + DABS(SR(LS))/2.0D0) IF (SMALL*ZTEST .NE. SMALL*TEST) GO TO 500 SR(LS) = 0.0D0 C ......EXIT GO TO 520 500 CONTINUE 510 CONTINUE 520 CONTINUE IF (LS .NE. L) GO TO 530 KASE = 3 GO TO 550 530 CONTINUE IF (LS .NE. M) GO TO 540 KASE = 1 GO TO 550 540 CONTINUE KASE = 2 L = LS 550 CONTINUE 560 CONTINUE L = L + 1 C C PERFORM THE TASK INDICATED BY KASE. C GO TO (570, 600, 620, 650), KASE C C DEFLATE NEGLIGIBLE SR(M). C 570 CONTINUE MM1 = M - 1 F = ER(M-1) ER(M-1) = 0.0D0 DO 590 KK = L, MM1 K = MM1 - KK + L T1 = SR(K) CALL RROTG(T1,F,CS,SN) SR(K) = T1 IF (K .EQ. L) GO TO 580 F = FLOP(-SN*ER(K-1)) ER(K-1) = FLOP(CS*ER(K-1)) 580 CONTINUE IF (WANTV) CALL RROT(P,VR(1,K),1,VR(1,M),1,CS,SN) IF (WANTV) CALL RROT(P,VI(1,K),1,VI(1,M),1,CS,SN) 590 CONTINUE GO TO 690 C C SPLIT AT NEGLIGIBLE SR(L). C 600 CONTINUE F = ER(L-1) ER(L-1) = 0.0D0 DO 610 K = L, M T1 = SR(K) CALL RROTG(T1,F,CS,SN) SR(K) = T1 F = FLOP(-SN*ER(K)) ER(K) = FLOP(CS*ER(K)) IF (WANTU) CALL RROT(N,UR(1,K),1,UR(1,L-1),1,CS,SN) IF (WANTU) CALL RROT(N,UI(1,K),1,UI(1,L-1),1,CS,SN) 610 CONTINUE GO TO 690 C C PERFORM ONE QR STEP. C 620 CONTINUE C C CALCULATE THE SHIFT. C SCALE = DMAX1(DABS(SR(M)),DABS(SR(M-1)),DABS(ER(M-1)), * DABS(SR(L)),DABS(ER(L))) SM = SR(M)/SCALE SMM1 = SR(M-1)/SCALE EMM1 = ER(M-1)/SCALE SL = SR(L)/SCALE EL = ER(L)/SCALE B = FLOP(((SMM1 + SM)*(SMM1 - SM) + EMM1**2)/2.0D0) C = FLOP((SM*EMM1)**2) SHIFT = 0.0D0 IF (B .EQ. 0.0D0 .AND. C .EQ. 0.0D0) GO TO 630 SHIFT = FLOP(DSQRT(B**2+C)) IF (B .LT. 0.0D0) SHIFT = -SHIFT SHIFT = FLOP(C/(B + SHIFT)) 630 CONTINUE F = FLOP((SL + SM)*(SL - SM) - SHIFT) G = FLOP(SL*EL) C C CHASE ZEROS. C MM1 = M - 1 DO 640 K = L, MM1 CALL RROTG(F,G,CS,SN) IF (K .NE. L) ER(K-1) = F F = FLOP(CS*SR(K) + SN*ER(K)) ER(K) = FLOP(CS*ER(K) - SN*SR(K)) G = FLOP(SN*SR(K+1)) SR(K+1) = FLOP(CS*SR(K+1)) IF (WANTV) CALL RROT(P,VR(1,K),1,VR(1,K+1),1,CS,SN) IF (WANTV) CALL RROT(P,VI(1,K),1,VI(1,K+1),1,CS,SN) CALL RROTG(F,G,CS,SN) SR(K) = F F = FLOP(CS*ER(K) + SN*SR(K+1)) SR(K+1) = FLOP(-SN*ER(K) + CS*SR(K+1)) G = FLOP(SN*ER(K+1)) ER(K+1) = FLOP(CS*ER(K+1)) IF (WANTU .AND. K .LT. N) * CALL RROT(N,UR(1,K),1,UR(1,K+1),1,CS,SN) IF (WANTU .AND. K .LT. N) * CALL RROT(N,UI(1,K),1,UI(1,K+1),1,CS,SN) 640 CONTINUE ER(M-1) = F ITER = ITER + 1 GO TO 690 C C CONVERGENCE C 650 CONTINUE C C MAKE THE SINGULAR VALUE POSITIVE C IF (SR(L) .GE. 0.0D0) GO TO 660 SR(L) = -SR(L) IF (WANTV) CALL WRSCAL(P,-1.0D0,VR(1,L),VI(1,L),1) 660 CONTINUE C C ORDER THE SINGULAR VALUE. C 670 IF (L .EQ. MM) GO TO 680 C ...EXIT IF (SR(L) .GE. SR(L+1)) GO TO 680 TR = SR(L) SR(L) = SR(L+1) SR(L+1) = TR IF (WANTV .AND. L .LT. P) * CALL WSWAP(P,VR(1,L),VI(1,L),1,VR(1,L+1),VI(1,L+1),1) IF (WANTU .AND. L .LT. N) * CALL WSWAP(N,UR(1,L),UI(1,L),1,UR(1,L+1),UI(1,L+1),1) L = L + 1 GO TO 670 680 CONTINUE ITER = 0 M = M - 1 690 CONTINUE GO TO 440 700 CONTINUE RETURN END SUBROUTINE WQRDC(XR,XI,LDX,N,P,QRAUXR,QRAUXI,JPVT,WORKR,WORKI, * JOB) INTEGER LDX,N,P,JOB INTEGER JPVT(1) DOUBLE PRECISION XR(LDX,1),XI(LDX,1),QRAUXR(1),QRAUXI(1), * WORKR(1),WORKI(1) C C WQRDC USES HOUSEHOLDER TRANSFORMATIONS TO COMPUTE THE QR C FACTORIZATION OF AN N BY P MATRIX X. COLUMN PIVOTING C BASED ON THE 2-NORMS OF THE REDUCED COLUMNS MAY BE C PERFORMED AT THE USERS OPTION. C C ON ENTRY C C X DOUBLE-COMPLEX(LDX,P), WHERE LDX .GE. N. C X CONTAINS THE MATRIX WHOSE DECOMPOSITION IS TO BE C COMPUTED. C C LDX INTEGER. C LDX IS THE LEADING DIMENSION OF THE ARRAY X. C C N INTEGER. C N IS THE NUMBER OF ROWS OF THE MATRIX X. C C P INTEGER. C P IS THE NUMBER OF COLUMNS OF THE MATRIX X. C C JPVT INTEGER(P). C JPVT CONTAINS INTEGERS THAT CONTROL THE SELECTION C OF THE PIVOT COLUMNS. THE K-TH COLUMN X(K) OF X C IS PLACED IN ONE OF THREE CLASSES ACCORDING TO THE C VALUE OF JPVT(K). C C IF JPVT(K) .GT. 0, THEN X(K) IS AN INITIAL C COLUMN. C C IF JPVT(K) .EQ. 0, THEN X(K) IS A FREE COLUMN. C C IF JPVT(K) .LT. 0, THEN X(K) IS A FINAL COLUMN. C C BEFORE THE DECOMPOSITION IS COMPUTED, INITIAL COLUMNS C ARE MOVED TO THE BEGINNING OF THE ARRAY X AND FINAL C COLUMNS TO THE END. BOTH INITIAL AND FINAL COLUMNS C ARE FROZEN IN PLACE DURING THE COMPUTATION AND ONLY C FREE COLUMNS ARE MOVED. AT THE K-TH STAGE OF THE C REDUCTION, IF X(K) IS OCCUPIED BY A FREE COLUMN C IT IS INTERCHANGED WITH THE FREE COLUMN OF LARGEST C REDUCED NORM. JPVT IS NOT REFERENCED IF C JOB .EQ. 0. C C WORK DOUBLE-COMPLEX(P). C WORK IS A WORK ARRAY. WORK IS NOT REFERENCED IF C JOB .EQ. 0. C C JOB INTEGER. C JOB IS AN INTEGER THAT INITIATES COLUMN PIVOTING. C IF JOB .EQ. 0, NO PIVOTING IS DONE. C IF JOB .NE. 0, PIVOTING IS DONE. C C ON RETURN C C X X CONTAINS IN ITS UPPER TRIANGLE THE UPPER C TRIANGULAR MATRIX R OF THE QR FACTORIZATION. C BELOW ITS DIAGONAL X CONTAINS INFORMATION FROM C WHICH THE UNITARY PART OF THE DECOMPOSITION C CAN BE RECOVERED. NOTE THAT IF PIVOTING HAS C BEEN REQUESTED, THE DECOMPOSITION IS NOT THAT C OF THE ORIGINAL MATRIX X BUT THAT OF X C WITH ITS COLUMNS PERMUTED AS DESCRIBED BY JPVT. C C QRAUX DOUBLE-COMPLEX(P). C QRAUX CONTAINS FURTHER INFORMATION REQUIRED TO RECOVER C THE UNITARY PART OF THE DECOMPOSITION. C C JPVT JPVT(K) CONTAINS THE INDEX OF THE COLUMN OF THE C ORIGINAL MATRIX THAT HAS BEEN INTERCHANGED INTO C THE K-TH COLUMN, IF PIVOTING WAS REQUESTED. C C LINPACK. THIS VERSION DATED 07/03/79 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C WQRDC USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. C C BLAS WAXPY,PYTHAG,WDOTCR,WDOTCI,WSCAL,WSWAP,WNRM2 C FORTRAN DABS,DIMAG,DMAX1,MIN0 C C INTERNAL VARIABLES C INTEGER J,JP,L,LP1,LUP,MAXJ,PL,PU DOUBLE PRECISION MAXNRM,WNRM2,TT DOUBLE PRECISION PYTHAG,WDOTCR,WDOTCI,NRMXLR,NRMXLI,TR,TI,FLOP LOGICAL NEGJ,SWAPJ C DOUBLE PRECISION ZDUMR,ZDUMI DOUBLE PRECISION CABS1 CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI) C PL = 1 PU = 0 IF (JOB .EQ. 0) GO TO 60 C C PIVOTING HAS BEEN REQUESTED. REARRANGE THE COLUMNS C ACCORDING TO JPVT. C DO 20 J = 1, P SWAPJ = JPVT(J) .GT. 0 NEGJ = JPVT(J) .LT. 0 JPVT(J) = J IF (NEGJ) JPVT(J) = -J IF (.NOT.SWAPJ) GO TO 10 IF (J .NE. PL) * CALL WSWAP(N,XR(1,PL),XI(1,PL),1,XR(1,J),XI(1,J),1) JPVT(J) = JPVT(PL) JPVT(PL) = J PL = PL + 1 10 CONTINUE 20 CONTINUE PU = P DO 50 JJ = 1, P J = P - JJ + 1 IF (JPVT(J) .GE. 0) GO TO 40 JPVT(J) = -JPVT(J) IF (J .EQ. PU) GO TO 30 CALL WSWAP(N,XR(1,PU),XI(1,PU),1,XR(1,J),XI(1,J),1) JP = JPVT(PU) JPVT(PU) = JPVT(J) JPVT(J) = JP 30 CONTINUE PU = PU - 1 40 CONTINUE 50 CONTINUE 60 CONTINUE C C COMPUTE THE NORMS OF THE FREE COLUMNS. C IF (PU .LT. PL) GO TO 80 DO 70 J = PL, PU QRAUXR(J) = WNRM2(N,XR(1,J),XI(1,J),1) QRAUXI(J) = 0.0D0 WORKR(J) = QRAUXR(J) WORKI(J) = QRAUXI(J) 70 CONTINUE 80 CONTINUE C C PERFORM THE HOUSEHOLDER REDUCTION OF X. C LUP = MIN0(N,P) DO 210 L = 1, LUP IF (L .LT. PL .OR. L .GE. PU) GO TO 120 C C LOCATE THE COLUMN OF LARGEST NORM AND BRING IT C INTO THE PIVOT POSITION. C MAXNRM = 0.0D0 MAXJ = L DO 100 J = L, PU IF (QRAUXR(J) .LE. MAXNRM) GO TO 90 MAXNRM = QRAUXR(J) MAXJ = J 90 CONTINUE 100 CONTINUE IF (MAXJ .EQ. L) GO TO 110 CALL WSWAP(N,XR(1,L),XI(1,L),1,XR(1,MAXJ),XI(1,MAXJ),1) QRAUXR(MAXJ) = QRAUXR(L) QRAUXI(MAXJ) = QRAUXI(L) WORKR(MAXJ) = WORKR(L) WORKI(MAXJ) = WORKI(L) JP = JPVT(MAXJ) JPVT(MAXJ) = JPVT(L) JPVT(L) = JP 110 CONTINUE 120 CONTINUE QRAUXR(L) = 0.0D0 QRAUXI(L) = 0.0D0 IF (L .EQ. N) GO TO 200 C C COMPUTE THE HOUSEHOLDER TRANSFORMATION FOR COLUMN L. C NRMXLR = WNRM2(N-L+1,XR(L,L),XI(L,L),1) NRMXLI = 0.0D0 IF (CABS1(NRMXLR,NRMXLI) .EQ. 0.0D0) GO TO 190 IF (CABS1(XR(L,L),XI(L,L)) .EQ. 0.0D0) GO TO 130 CALL WSIGN(NRMXLR,NRMXLI,XR(L,L),XI(L,L),NRMXLR,NRMXLI) 130 CONTINUE CALL WDIV(1.0D0,0.0D0,NRMXLR,NRMXLI,TR,TI) CALL WSCAL(N-L+1,TR,TI,XR(L,L),XI(L,L),1) XR(L,L) = FLOP(1.0D0 + XR(L,L)) C C APPLY THE TRANSFORMATION TO THE REMAINING COLUMNS, C UPDATING THE NORMS. C LP1 = L + 1 IF (P .LT. LP1) GO TO 180 DO 170 J = LP1, P TR = -WDOTCR(N-L+1,XR(L,L),XI(L,L),1,XR(L,J), * XI(L,J),1) TI = -WDOTCI(N-L+1,XR(L,L),XI(L,L),1,XR(L,J), * XI(L,J),1) CALL WDIV(TR,TI,XR(L,L),XI(L,L),TR,TI) CALL WAXPY(N-L+1,TR,TI,XR(L,L),XI(L,L),1,XR(L,J), * XI(L,J),1) IF (J .LT. PL .OR. J .GT. PU) GO TO 160 IF (CABS1(QRAUXR(J),QRAUXI(J)) .EQ. 0.0D0) * GO TO 160 TT = 1.0D0 - (PYTHAG(XR(L,J),XI(L,J))/QRAUXR(J))**2 TT = DMAX1(TT,0.0D0) TR = FLOP(TT) TT = FLOP(1.0D0+0.05D0*TT*(QRAUXR(J)/WORKR(J))**2) IF (TT .EQ. 1.0D0) GO TO 140 QRAUXR(J) = QRAUXR(J)*DSQRT(TR) QRAUXI(J) = QRAUXI(J)*DSQRT(TR) GO TO 150 140 CONTINUE QRAUXR(J) = WNRM2(N-L,XR(L+1,J),XI(L+1,J),1) QRAUXI(J) = 0.0D0 WORKR(J) = QRAUXR(J) WORKI(J) = QRAUXI(J) 150 CONTINUE 160 CONTINUE 170 CONTINUE 180 CONTINUE C C SAVE THE TRANSFORMATION. C QRAUXR(L) = XR(L,L) QRAUXI(L) = XI(L,L) XR(L,L) = -NRMXLR XI(L,L) = -NRMXLI 190 CONTINUE 200 CONTINUE 210 CONTINUE RETURN END SUBROUTINE WQRSL(XR,XI,LDX,N,K,QRAUXR,QRAUXI,YR,YI,QYR,QYI,QTYR, * QTYI,BR,BI,RSDR,RSDI,XBR,XBI,JOB,INFO) INTEGER LDX,N,K,JOB,INFO DOUBLE PRECISION XR(LDX,1),XI(LDX,1),QRAUXR(1),QRAUXI(1),YR(1), * YI(1),QYR(1),QYI(1),QTYR(1),QTYI(1),BR(1),BI(1), * RSDR(1),RSDI(1),XBR(1),XBI(1) C C WQRSL APPLIES THE OUTPUT OF WQRDC TO COMPUTE COORDINATE C TRANSFORMATIONS, PROJECTIONS, AND LEAST SQUARES SOLUTIONS. C FOR K .LE. MIN(N,P), LET XK BE THE MATRIX C C XK = (X(JPVT(1)),X(JPVT(2)), ... ,X(JPVT(K))) C C FORMED FROM COLUMNNS JPVT(1), ... ,JPVT(K) OF THE ORIGINAL C N X P MATRIX X THAT WAS INPUT TO WQRDC (IF NO PIVOTING WAS C DONE, XK CONSISTS OF THE FIRST K COLUMNS OF X IN THEIR C ORIGINAL ORDER). WQRDC PRODUCES A FACTORED UNITARY MATRIX Q C AND AN UPPER TRIANGULAR MATRIX R SUCH THAT C C XK = Q * (R) C (0) C C THIS INFORMATION IS CONTAINED IN CODED FORM IN THE ARRAYS C X AND QRAUX. C C ON ENTRY C C X DOUBLE-COMPLEX(LDX,P). C X CONTAINS THE OUTPUT OF WQRDC. C C LDX INTEGER. C LDX IS THE LEADING DIMENSION OF THE ARRAY X. C C N INTEGER. C N IS THE NUMBER OF ROWS OF THE MATRIX XK. IT MUST C HAVE THE SAME VALUE AS N IN WQRDC. C C K INTEGER. C K IS THE NUMBER OF COLUMNS OF THE MATRIX XK. K C MUST NNOT BE GREATER THAN MIN(N,P), WHERE P IS THE C SAME AS IN THE CALLING SEQUENCE TO WQRDC. C C QRAUX DOUBLE-COMPLEX(P). C QRAUX CONTAINS THE AUXILIARY OUTPUT FROM WQRDC. C C Y DOUBLE-COMPLEX(N) C Y CONTAINS AN N-VECTOR THAT IS TO BE MANIPULATED C BY WQRSL. C C JOB INTEGER. C JOB SPECIFIES WHAT IS TO BE COMPUTED. JOB HAS C THE DECIMAL EXPANSION ABCDE, WITH THE FOLLOWING C MEANING. C C IF A.NE.0, COMPUTE QY. C IF B,C,D, OR E .NE. 0, COMPUTE QTY. C IF C.NE.0, COMPUTE B. C IF D.NE.0, COMPUTE RSD. C IF E.NE.0, COMPUTE XB. C C NOTE THAT A REQUEST TO COMPUTE B, RSD, OR XB C AUTOMATICALLY TRIGGERS THE COMPUTATION OF QTY, FOR C WHICH AN ARRAY MUST BE PROVIDED IN THE CALLING C SEQUENCE. C C ON RETURN C C QY DOUBLE-COMPLEX(N). C QY CONNTAINS Q*Y, IF ITS COMPUTATION HAS BEEN C REQUESTED. C C QTY DOUBLE-COMPLEX(N). C QTY CONTAINS CTRANS(Q)*Y, IF ITS COMPUTATION HAS C BEEN REQUESTED. HERE CTRANS(Q) IS THE CONJUGATE C TRANSPOSE OF THE MATRIX Q. C C B DOUBLE-COMPLEX(K) C B CONTAINS THE SOLUTION OF THE LEAST SQUARES PROBLEM C C MINIMIZE NORM2(Y - XK*B), C C IF ITS COMPUTATION HAS BEEN REQUESTED. (NOTE THAT C IF PIVOTING WAS REQUESTED IN WQRDC, THE J-TH C COMPONENT OF B WILL BE ASSOCIATED WITH COLUMN JPVT(J) C OF THE ORIGINAL MATRIX X THAT WAS INPUT INTO WQRDC.) C C RSD DOUBLE-COMPLEX(N). C RSD CONTAINS THE LEAST SQUARES RESIDUAL Y - XK*B, C IF ITS COMPUTATION HAS BEEN REQUESTED. RSD IS C ALSO THE ORTHOGONAL PROJECTION OF Y ONTO THE C ORTHOGONAL COMPLEMENT OF THE COLUMN SPACE OF XK. C C XB DOUBLE-COMPLEX(N). C XB CONTAINS THE LEAST SQUARES APPROXIMATION XK*B, C IF ITS COMPUTATION HAS BEEN REQUESTED. XB IS ALSO C THE ORTHOGONAL PROJECTION OF Y ONTO THE COLUMN SPACE C OF X. C C INFO INTEGER. C INFO IS ZERO UNLESS THE COMPUTATION OF B HAS C BEEN REQUESTED AND R IS EXACTLY SINGULAR. IN C THIS CASE, INFO IS THE INDEX OF THE FIRST ZERO C DIAGONAL ELEMENT OF R AND B IS LEFT UNALTERED. C C THE PARAMETERS QY, QTY, B, RSD, AND XB ARE NOT REFERENCED C IF THEIR COMPUTATION IS NOT REQUESTED AND IN THIS CASE C CAN BE REPLACED BY DUMMY VARIABLES IN THE CALLING PROGRAM. C TO SAVE STORAGE, THE USER MAY IN SOME CASES USE THE SAME C ARRAY FOR DIFFERENT PARAMETERS IN THE CALLING SEQUENCE. A C FREQUENTLY OCCURING EXAMPLE IS WHEN ONE WISHES TO COMPUTE C ANY OF B, RSD, OR XB AND DOES NOT NEED Y OR QTY. IN THIS C CASE ONE MAY IDENTIFY Y, QTY, AND ONE OF B, RSD, OR XB, WHILE C PROVIDING SEPARATE ARRAYS FOR ANYTHING ELSE THAT IS TO BE C COMPUTED. THUS THE CALLING SEQUENCE C C CALL WQRSL(X,LDX,N,K,QRAUX,Y,DUM,Y,B,Y,DUM,110,INFO) C C WILL RESULT IN THE COMPUTATION OF B AND RSD, WITH RSD C OVERWRITING Y. MORE GENERALLY, EACH ITEM IN THE FOLLOWING C LIST CONTAINS GROUPS OF PERMISSIBLE IDENTIFICATIONS FOR C A SINGLE CALLINNG SEQUENCE. C C 1. (Y,QTY,B) (RSD) (XB) (QY) C C 2. (Y,QTY,RSD) (B) (XB) (QY) C C 3. (Y,QTY,XB) (B) (RSD) (QY) C C 4. (Y,QY) (QTY,B) (RSD) (XB) C C 5. (Y,QY) (QTY,RSD) (B) (XB) C C 6. (Y,QY) (QTY,XB) (B) (RSD) C C IN ANY GROUP THE VALUE RETURNED IN THE ARRAY ALLOCATED TO C THE GROUP CORRESPONDS TO THE LAST MEMBER OF THE GROUP. C C LINPACK. THIS VERSION DATED 07/03/79 . C G.W. STEWART, UNIVERSITY OF MARYLAND, ARGONNE NATIONAL LAB. C C WQRSL USES THE FOLLOWING FUNCTIONS AND SUBPROGRAMS. C C BLAS WAXPY,WCOPY,WDOTCR,WDOTCI C FORTRAN DABS,DIMAG,MIN0,MOD C C INTERNAL VARIABLES C INTEGER I,J,JJ,JU,KP1 DOUBLE PRECISION WDOTCR,WDOTCI,TR,TI,TEMPR,TEMPI LOGICAL CB,CQY,CQTY,CR,CXB C DOUBLE PRECISION ZDUMR,ZDUMI DOUBLE PRECISION CABS1 CABS1(ZDUMR,ZDUMI) = DABS(ZDUMR) + DABS(ZDUMI) C C SET INFO FLAG. C INFO = 0 C C DETERMINE WHAT IS TO BE COMPUTED. C CQY = JOB/10000 .NE. 0 CQTY = MOD(JOB,10000) .NE. 0 CB = MOD(JOB,1000)/100 .NE. 0 CR = MOD(JOB,100)/10 .NE. 0 CXB = MOD(JOB,10) .NE. 0 JU = MIN0(K,N-1) C C SPECIAL ACTION WHEN N=1. C IF (JU .NE. 0) GO TO 80 IF (.NOT.CQY) GO TO 10 QYR(1) = YR(1) QYI(1) = YI(1) 10 CONTINUE IF (.NOT.CQTY) GO TO 20 QTYR(1) = YR(1) QTYI(1) = YI(1) 20 CONTINUE IF (.NOT.CXB) GO TO 30 XBR(1) = YR(1) XBI(1) = YI(1) 30 CONTINUE IF (.NOT.CB) GO TO 60 IF (CABS1(XR(1,1),XI(1,1)) .NE. 0.0D0) GO TO 40 INFO = 1 GO TO 50 40 CONTINUE CALL WDIV(YR(1),YI(1),XR(1,1),XI(1,1),BR(1),BI(1)) 50 CONTINUE 60 CONTINUE IF (.NOT.CR) GO TO 70 RSDR(1) = 0.0D0 RSDI(1) = 0.0D0 70 CONTINUE GO TO 290 80 CONTINUE C C SET UP TO COMPUTE QY OR QTY. C IF (CQY) CALL WCOPY(N,YR,YI,1,QYR,QYI,1) IF (CQTY) CALL WCOPY(N,YR,YI,1,QTYR,QTYI,1) IF (.NOT.CQY) GO TO 110 C C COMPUTE QY. C DO 100 JJ = 1, JU J = JU - JJ + 1 IF (CABS1(QRAUXR(J),QRAUXI(J)) .EQ. 0.0D0) * GO TO 90 TEMPR = XR(J,J) TEMPI = XI(J,J) XR(J,J) = QRAUXR(J) XI(J,J) = QRAUXI(J) TR = -WDOTCR(N-J+1,XR(J,J),XI(J,J),1,QYR(J),QYI(J),1) TI = -WDOTCI(N-J+1,XR(J,J),XI(J,J),1,QYR(J),QYI(J),1) CALL WDIV(TR,TI,XR(J,J),XI(J,J),TR,TI) CALL WAXPY(N-J+1,TR,TI,XR(J,J),XI(J,J),1,QYR(J), * QYI(J),1) XR(J,J) = TEMPR XI(J,J) = TEMPI 90 CONTINUE 100 CONTINUE 110 CONTINUE IF (.NOT.CQTY) GO TO 140 C C COMPUTE CTRANS(Q)*Y. C DO 130 J = 1, JU IF (CABS1(QRAUXR(J),QRAUXI(J)) .EQ. 0.0D0) * GO TO 120 TEMPR = XR(J,J) TEMPI = XI(J,J) XR(J,J) = QRAUXR(J) XI(J,J) = QRAUXI(J) TR = -WDOTCR(N-J+1,XR(J,J),XI(J,J),1,QTYR(J), * QTYI(J),1) TI = -WDOTCI(N-J+1,XR(J,J),XI(J,J),1,QTYR(J), * QTYI(J),1) CALL WDIV(TR,TI,XR(J,J),XI(J,J),TR,TI) CALL WAXPY(N-J+1,TR,TI,XR(J,J),XI(J,J),1,QTYR(J), * QTYI(J),1) XR(J,J) = TEMPR XI(J,J) = TEMPI 120 CONTINUE 130 CONTINUE 140 CONTINUE C C SET UP TO COMPUTE B, RSD, OR XB. C IF (CB) CALL WCOPY(K,QTYR,QTYI,1,BR,BI,1) KP1 = K + 1 IF (CXB) CALL WCOPY(K,QTYR,QTYI,1,XBR,XBI,1) IF (CR .AND. K .LT. N) * CALL WCOPY(N-K,QTYR(KP1),QTYI(KP1),1,RSDR(KP1),RSDI(KP1),1) IF (.NOT.CXB .OR. KP1 .GT. N) GO TO 160 DO 150 I = KP1, N XBR(I) = 0.0D0 XBI(I) = 0.0D0 150 CONTINUE 160 CONTINUE IF (.NOT.CR) GO TO 180 DO 170 I = 1, K RSDR(I) = 0.0D0 RSDI(I) = 0.0D0 170 CONTINUE 180 CONTINUE IF (.NOT.CB) GO TO 230 C C COMPUTE B. C DO 210 JJ = 1, K J = K - JJ + 1 IF (CABS1(XR(J,J),XI(J,J)) .NE. 0.0D0) GO TO 190 INFO = J C ......EXIT C ......EXIT GO TO 220 190 CONTINUE CALL WDIV(BR(J),BI(J),XR(J,J),XI(J,J),BR(J),BI(J)) IF (J .EQ. 1) GO TO 200 TR = -BR(J) TI = -BI(J) CALL WAXPY(J-1,TR,TI,XR(1,J),XI(1,J),1,BR,BI,1) 200 CONTINUE 210 CONTINUE 220 CONTINUE 230 CONTINUE IF (.NOT.CR .AND. .NOT.CXB) GO TO 280 C C COMPUTE RSD OR XB AS REQUIRED. C DO 270 JJ = 1, JU J = JU - JJ + 1 IF (CABS1(QRAUXR(J),QRAUXI(J)) .EQ. 0.0D0) * GO TO 260 TEMPR = XR(J,J) TEMPI = XI(J,J) XR(J,J) = QRAUXR(J) XI(J,J) = QRAUXI(J) IF (.NOT.CR) GO TO 240 TR = -WDOTCR(N-J+1,XR(J,J),XI(J,J),1,RSDR(J), * RSDI(J),1) TI = -WDOTCI(N-J+1,XR(J,J),XI(J,J),1,RSDR(J), * RSDI(J),1) CALL WDIV(TR,TI,XR(J,J),XI(J,J),TR,TI) CALL WAXPY(N-J+1,TR,TI,XR(J,J),XI(J,J),1,RSDR(J), * RSDI(J),1) 240 CONTINUE IF (.NOT.CXB) GO TO 250 TR = -WDOTCR(N-J+1,XR(J,J),XI(J,J),1,XBR(J), * XBI(J),1) TI = -WDOTCI(N-J+1,XR(J,J),XI(J,J),1,XBR(J), * XBI(J),1) CALL WDIV(TR,TI,XR(J,J),XI(J,J),TR,TI) CALL WAXPY(N-J+1,TR,TI,XR(J,J),XI(J,J),1,XBR(J), * XBI(J),1) 250 CONTINUE XR(J,J) = TEMPR XI(J,J) = TEMPI 260 CONTINUE 270 CONTINUE 280 CONTINUE 290 CONTINUE RETURN END SUBROUTINE MAGIC(A,LDA,N) C C ALGORITHMS FOR MAGIC SQUARES TAKEN FROM C MATHEMATICAL RECREATIONS AND ESSAYS, 12TH ED., C BY W. W. ROUSE BALL AND H. S. M. COXETER C DOUBLE PRECISION A(LDA,N),T C IF (MOD(N,4) .EQ. 0) GO TO 100 IF (MOD(N,2) .EQ. 0) M = N/2 IF (MOD(N,2) .NE. 0) M = N C C ODD ORDER OR UPPER CORNER OF EVEN ORDER C DO 20 J = 1,M DO 10 I = 1,M A(I,J) = 0 10 CONTINUE 20 CONTINUE I = 1 J = (M+1)/2 MM = M*M DO 40 K = 1, MM A(I,J) = K I1 = I-1 J1 = J+1 IF(I1.LT.1) I1 = M IF(J1.GT.M) J1 = 1 IF(IDINT(A(I1,J1)).EQ.0) GO TO 30 I1 = I+1 J1 = J 30 I = I1 J = J1 40 CONTINUE IF (MOD(N,2) .NE. 0) RETURN C C REST OF EVEN ORDER C T = M*M DO 60 I = 1, M DO 50 J = 1, M IM = I+M JM = J+M A(I,JM) = A(I,J) + 2*T A(IM,J) = A(I,J) + 3*T A(IM,JM) = A(I,J) + T 50 CONTINUE 60 CONTINUE M1 = (M-1)/2 IF (M1.EQ.0) RETURN DO 70 J = 1, M1 CALL RSWAP(M,A(1,J),1,A(M+1,J),1) 70 CONTINUE M1 = (M+1)/2 M2 = M1 + M CALL RSWAP(1,A(M1,1),1,A(M2,1),1) CALL RSWAP(1,A(M1,M1),1,A(M2,M1),1) M1 = N+1-(M-3)/2 IF(M1.GT.N) RETURN DO 80 J = M1, N CALL RSWAP(M,A(1,J),1,A(M+1,J),1) 80 CONTINUE RETURN C C DOUBLE EVEN ORDER C 100 K = 1 DO 120 I = 1, N DO 110 J = 1, N A(I,J) = K IF (MOD(I,4)/2 .EQ. MOD(J,4)/2) A(I,J) = N*N+1 - K K = K+1 110 CONTINUE 120 CONTINUE RETURN END SUBROUTINE BASE(X,B,EPS,S,N) DOUBLE PRECISION X,B,EPS,S(1),T C C STORE BASE B REPRESENTATION OF X IN S(1:N) C INTEGER PLUS,MINUS,DOT,ZERO,COMMA DATA PLUS/41/,MINUS/42/,DOT/47/,ZERO/0/,COMMA/48/ L = 1 IF (X .GE. 0.0D0) S(L) = PLUS IF (X .LT. 0.0D0) S(L) = MINUS S(L+1) = ZERO S(L+2) = DOT X = DABS(X) IF (X .NE. 0.0D0) K = DLOG(X)/DLOG(B) IF (X .EQ. 0.0D0) K = 0 IF (X .GT. 1.0D0) K = K + 1 X = X/B**K IF (B*X .GE. B) K = K + 1 IF (B*X .GE. B) X = X/B IF (EPS .NE. 0.0D0) M = -DLOG(EPS)/DLOG(B) + 4 IF (EPS .EQ. 0.0D0) M = 54 DO 10 L = 4, M X = B*X J = IDINT(X) S(L) = DFLOAT(J) X = X - S(L) 10 CONTINUE S(M+1) = COMMA IF (K .GE. 0) S(M+2) = PLUS IF (K .LT. 0) S(M+2) = MINUS T = DABS(DFLOAT(K)) N = M + 3 IF (T .GE. B) N = N + IDINT(DLOG(T)/DLOG(B)) L = N 20 J = IDINT(DMOD(T,B)) S(L) = DFLOAT(J) L = L - 1 T = T/B IF (L .GE. M+3) GO TO 20 RETURN END DOUBLE PRECISION FUNCTION URAND(IY) INTEGER IY C C URAND IS A UNIFORM RANDOM NUMBER GENERATOR BASED ON THEORY AND C SUGGESTIONS GIVEN IN D.E. KNUTH (1969), VOL 2. THE INTEGER IY C SHOULD BE INITIALIZED TO AN ARBITRARY INTEGER PRIOR TO THE FIRST CALL C TO URAND. THE CALLING PROGRAM SHOULD NOT ALTER THE VALUE OF IY C BETWEEN SUBSEQUENT CALLS TO URAND. VALUES OF URAND WILL BE RETURNED C IN THE INTERVAL (0,1). C INTEGER IA,IC,ITWO,M2,M,MIC DOUBLE PRECISION HALFM,S DOUBLE PRECISION DATAN,DSQRT DATA M2/0/,ITWO/2/ IF (M2 .NE. 0) GO TO 20 C C IF FIRST ENTRY, COMPUTE MACHINE INTEGER WORD LENGTH C M = 1 10 M2 = M M = ITWO*M2 IF (M .GT. M2) GO TO 10 HALFM = M2 C C COMPUTE MULTIPLIER AND INCREMENT FOR LINEAR CONGRUENTIAL METHOD C IA = 8*IDINT(HALFM*DATAN(1.D0)/8.D0) + 5 IC = 2*IDINT(HALFM*(0.5D0-DSQRT(3.D0)/6.D0)) + 1 MIC = (M2 - IC) + M2 C C S IS THE SCALE FACTOR FOR CONVERTING TO FLOATING POINT C S = 0.5D0/HALFM C C COMPUTE NEXT RANDOM NUMBER C 20 IY = IY*IA C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHICH DO NOT ALLOW C INTEGER OVERFLOW ON ADDITION C IF (IY .GT. MIC) IY = (IY - M2) - M2 C IY = IY + IC C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE THE C WORD LENGTH FOR ADDITION IS GREATER THAN FOR MULTIPLICATION C IF (IY/2 .GT. M2) IY = (IY - M2) - M2 C C THE FOLLOWING STATEMENT IS FOR COMPUTERS WHERE INTEGER C OVERFLOW AFFECTS THE SIGN BIT C IF (IY .LT. 0) IY = (IY + M2) + M2 URAND = DFLOAT(IY)*S RETURN END SUBROUTINE WMUL(AR,AI,BR,BI,CR,CI) DOUBLE PRECISION AR,AI,BR,BI,CR,CI,T,FLOP C C = A*B T = AR*BI + AI*BR IF (T .NE. 0.0D0) T = FLOP(T) CR = FLOP(AR*BR - AI*BI) CI = T RETURN END SUBROUTINE WDIV(AR,AI,BR,BI,CR,CI) DOUBLE PRECISION AR,AI,BR,BI,CR,CI C C = A/B DOUBLE PRECISION S,D,ARS,AIS,BRS,BIS,FLOP S = DABS(BR) + DABS(BI) IF (S .EQ. 0.0D0) CALL ERROR(27) IF (S .EQ. 0.0D0) RETURN ARS = AR/S AIS = AI/S BRS = BR/S BIS = BI/S D = BRS**2 + BIS**2 SHAR_EOF # End of shell archive exit 0 -- Bob Page, U of Lowell CS Dept. page@swan.ulowell.edu ulowell!page Have five nice days.